# -*- coding: utf-8 -*-

import numpy as np;
import matplotlib.pyplot as mp;

def step_euler(y,t,h,f):
	return y + h * f(y, t);

def step_heun(y, t, h, f):
    pn1 = f(y, t)
    yn1 = y + h * pn1
    pn2 = f(yn1, t + h)
    return y + ((h * (pn1 + pn2)) / 2)
    

def step_middle(y,t,h,f):
	y05 = y + h/2 * f(y,t);
	pn = f(y05, t+h/2);
	return y + h * pn;


def step_RK4(y,t,h,f):
        pn1 = f(y, t)
        yn1 = y + (1/2) * h * pn1
        pn2 = f(yn1, t + (1/2) * h)
        yn2 = y + (1/2) * h * pn2
        pn3 = f(yn2, t + (1/2)*h)
        yn3 = y + h * pn3
        pn4 = f(yn3, t + h)
        
	return y + ((h * (pn1 + 2*pn2 + 2*pn3 + pn4))/6);
	

def meth_n_step(y0,t0,N,h,f,meth):
	y = y0;
	t = t0;
	
	for i in np.arange(0, N):
		y = meth(y, t, h, f);
		t = t + h;
	return y;
	
def meth_epsilon(y0,t0,tf,eps,f,meth):
	N = 1;
	h = (tf - t0) / N;
	y = meth_n_step(y0,t0,N,h,f,meth);
	
	while(True):
		N = N * 2;
		h = (tf - t0) / N;
		yl = y;
		y = meth_n_step(y0,t0,N,h,f,meth);
		if (np.linalg.norm(y - yl) < eps):
			print "Precision atteinte pour "+f.__name__+". Il a fallut "+str(N)+" pas";
			break;
	return y;
	
def compute_step(y0,t0,tf,eps,f,meth):
	N = 1;
	h = (tf - t0) / N;
	y = meth_n_step(y0,t0,N,h,f,meth);
	
	while(True):
		N = N * 2;
		h = (tf - t0) / N;
		yl = y;
		y = meth_n_step(y0,t0,N,h,f,meth);
		if (np.linalg.norm(y - yl) < eps):
			#~ print "Precision atteinte pour "+f.__name__+". Il a fallut "+str(N)+" pas";
			break;
	return N;
	
def fill_graph(cauchy, X, eps, meth):
	n = len(X);
	if N<n:
		N=n;
	#~ print str(N) + " pas sont necessaires.";
	
	h = (X[n-1] - cauchy[1]) / N;
	Y = np.zeros([len(X),cauchy[0]]);
		
	y = cauchy[2];
	t = cauchy[1];
	x = X[0];
	
	j = 0;
	i = 0;
	while j < len(X):
		if t>=X[j]:
			Y[j] = y;
			j+=1;
		y = meth(y, t, h, cauchy[3]);
		yl = meth(y, t, h/2, cauchy[3]);
		yl = meth(y, t+h/2, h/2, cauchy[3]);
		while (np.linalg.norm(y - yl) > eps and N < 200000): # Augmentation du pas si precision insuffisante. (limitation du pas a 10000)
			N = N*2;
			print "augmentation a "+str(N)
			h = h/2;
			yl = meth(y, t, h, cauchy[3]);
		t = t + h;
		i+=1;
		
	return Y;
	
def fill_graph_noopt(cauchy, X, eps, meth):
	n = len(X);
	N = compute_step(cauchy[2], cauchy[1], X[n-1], eps, cauchy[3], meth);
	if N<n:
		N=n;
	#~ print str(N) + " pas sont necessaires.";
	
	h = (X[n-1] - cauchy[1]) / N;
	Y = np.zeros([len(X),cauchy[0]]);
		
	y = cauchy[2];
	t = cauchy[1];
	x = X[0];
	
	j = 0;
	i = 0;
	while j < len(X):
		if t>=X[j]:
			Y[j] = y;
			j+=1;
		y = meth(y, t, h, cauchy[3]);
		t = t + h;
		print t;
		i+=1;
		
	return Y;
	
	
# Cette fonction essaie de calculer la periode d'une fonction.
# Elle calcule la moyenne et compte le nombre de traversées montantes
def compute_period(X,Y):
	if len(X) != len(Y):
		return -1;
	n = len(Y);
	mean = sum(Y) / n;
	nup = 0;
	
	
	for i in np.arange(0, n-1):
		if Y[i]<=mean and mean<=Y[i+1]:
			nup +=1;
			
	if nup == 0:
		print "Pas de periodicité";
		return -1;
	return (X[len(X)-1] - X[0]) / nup;
		
		
def trace_local(cauchy, diff):
	mp.clf();
	X = np.arange(0,1,0.01);
	variation = np.arange(-diff,diff,diff *2/8);
	for var in variation:
		tmp_cauchy = (cauchy[0], cauchy[1], cauchy[2]+var, cauchy[3])
		Y = fill_graph(tmp_cauchy, X, 0.05, step_euler);
		mp.plot(X, Y[:,0]);
	mp.title("Solutions proches de y0 = "+str(cauchy[2]));
